function p = spmitrf(m)
% spmitrf  - uncompile an affine transformation into paramters
%
% FORMAT:       p = spmitrf(m)
%
% Input fields:
%
%       m           4x4 affine transformation matrix
%
% Output fields:
%
%       p           1x4 cell array for use with spmtrf

% Version:  v0.7g
% Build:    9031721
% Date:     Mar-17 2009, 9:37 PM CET
% Author:   John Ashburner & Stefan Kiebel, Wellcome Department of I.N.
% Editor:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin ~= 1 || ...
   ~isa(m, 'double') || ...
   ~isequal(size(m), [4, 4]) || ...
    any(isinf(m(:)) | isnan(m(:))) || ...
    any(m(4, :) ~= [0, 0, 0, 1])
    error( ...
        'BVQXtools:BadArgument', ...
        'Invalid or missing transformation matrix supplied.' ...
    );
end
p = cell(1, 4);

% translation
p{1} = m(1:3, 4)';

% decompose rotation matrix into zooms
p{2} = [0, 0, 0];
r = m(1:3, 1:3);
rc = chol(r' * r);
p{3} = diag(rc)';
if det(rc) < 0
    p{3}(1) = -p{3}(1);
end

% shear
rc = diag(diag(rc)) \ rc;
p{4} = rc([4, 7, 8]);

% create untranslated for rotation
r0 = spmtrf([0, 0, 0], p{2:4});
r0 = r0(1:3, 1:3);
r0 = r / r0;

% compute rotation parameters
p{2}(2)  = asin(rang(r0(1, 3)));
if (abs(p{2}(2)) - pi / 2) ^ 2 < 1e-9
	p{2}(1) = 0;
	p{2}(3) = atan2(-rang(r0(2, 1)), rang(-r0(3, 1) / r0(1, 3)));
else
	c = cos(p{2}(2));
	p{2}(1) = atan2(rang(r0(2, 3) / c), rang(r0(3, 3) / c));
	p{2}(3) = atan2(rang(r0(1, 2) / c), rang(r0(1, 1) / c));
end

% There may be slight rounding errors making b>1 or b<-1.
function a = rang(b)
a = min(max(b, -1), 1);
